;+
; NAME:
;       MAKE_RADFIELD_PLOT_PAPER
;
; PURPOSE:
;       Make the radiation field comparison plots for the IRDC paper
;
; CATEGORY:
;       Paper plot creation
;
; OPTIONAL INPUTS:
;       SUBSET -- Subset to use
;
;-

PRO MAKE_RADFIELD_PLOT_PAPER, SUBSET=subset
  
  IF n_elements(subset) EQ 0 THEN subset = ''
  
  myps,'./emaf_paper/plots/plot_radfield_paper.eps',/bw,xsize=3.1,ysize=3.1
  cs = 0.7
  xrange = [2,12]
  leg_cs = 0.85*cs
  GOTO,bottompanel
  
  ;; Define the COMMON BLOCK
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  ;; Since we are using routines buried in BGPS_FFORE, compile now to
  ;;   be sure they are available.
  ;; FORWARD_FUNCTION BGPS_FFORE_MODEL, BGPS_FFORE_MODEL2, BGPS_FFORE_SIMPLE, $
  ;;    BGPS_FFORE_SIMPLE2, BGPS_FFORE_LOOPIE, STAR_DENSITY
  ;; RESOLVE_ROUTINE, 'BGPS_FFORE', /COMPILE_FULL_FILE, /EITHER, /NO_RECOMPILE
  
  IF n_elements(rb3) EQ 0 THEN $
     restore,'./irdc_dist_model/bgps_rb3'+subset+'.sav',/ver
  n = n_elements(rb3)
  
  ;; Define arrays for TEMP_DEPENDENT_VARIABLES
  tau    = dblarr(n)
  f_data = dblarr(n)
  sig_f  = dblarr(n)
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  
  ;; Using !MW parameters -- Goes into COMMON block!
  R0 = !MW.R0
  d = dindgen(!MW.NBINS*5.)*!MW.BINSIZE + !MW.BINSTART
  corr = R0 / 8.5d3
  
  tstring = 'FIT' 
  
  ;; There are currently 3 cases of the model
  nmods     = 3
  dfn       = strarr(nmods)
  thisnf    = bytarr(n)
  fitdist   = fltarr(n)
  ffore_vec = fltarr(n)
  
  ;; Set up structures for the radiation fields
  radfield = replicate( create_struct('rgal', fltarr(!MW.NBINS*5.), $
                                      'z', fltarr(!MW.NBINS*5.)), nmods)
  
  myps,'./emaf_paper/plots/plot_radfield_paper.eps',/bw,xsize=3.1,ysize=3.1
  
  cs = 0.7
  xrange = [2,12]
  leg_cs = 0.85*cs
  GOTO,bottompanel
  ;;=========================================================================
  ;; TOP PANEL -- Cases of the radiation field in the galaxy
  
  multiplot_xm,/xmargin,[1,2]
  leg_cs = 0.85*cs
  srho = cgSymbol('rho')
  
  ;; cnames = ['','','','YGB7','Lime Green','TG3','RED7']
  cnames = ['YGB6','Green','TG3']
  colors = cgColor(cnames)
  ;; textcolor = cgColor('RYB8')
  textcolor = cgColor('Opposite')
  print,colors
  thick=[6,5,4]
  
  do_Tfit = 1b
  
  rotc=0b
  
  ;; Do re-usable calculations
  GENERATE_LPSTR, ROTC=rotc, SUFF=suff
  
  xtit = 'R!dgal!n [kpc]  (Z = 0)'
  ;; ytit = cgSymbol('rho')+'(R!dgal!n) / '+cgSymbol('rho')+'(R=R!d0!n)' 
  ytit = 'Normalized Model Volume Density at Z=0'
  ytf = 'exponent10'            ;"(F0.1)"
  
  yst=9
  
  plot,[0,1],[0,1],ytickformat=ytf,charsize=1.2,xr=xrange,$
       /xst,yr=[1.d-1,2.d1],yst=yst,/NODATA,ytit=ytit,$;xticks=3,xminor=5,$
       xmargin=[10,10],/ylog,color=cgColor('Opposite')
  
  chmin = fltarr(nmods)
  ndof  = fltarr(nmods)
  tdp   = fltarr(nmods)
  
  FOR j=0, nmods-1 DO BEGIN
     
     cj = j
     
     message,'Case #'+string(j+1,format="(I0)"),/inf
     dfn[j] = './irdc_dist_model/data/case'+$
              string(j+1,format="(I0)")+$
              '_irdc_TdFIT'+subset+suff+'.sav'
     
     restore,dfn[j]
     chmin[j] = bestnorm
     ndof[j] = dof
     
     message,'Minimum Chisq = '+string(bestnorm,format="(F0.2)"),/inf
     
     ;; First, do a vector along Z=0, from the GC outward & normalize
     ;; Get the normalization (rhon)
     R = [R0]
     Z = [0.d]
     rho_hi = hi_density(R,Z)
     rho_h2 = h2_density(R,Z)
     rho_star = star_density(R,Z)
     rhin = rho_hi
     rh2n = rho_h2
     rstn = rho_star
     
     CASE res[5] OF
        0: fj = BGPS_FFORE_MODEL(R0,res[0],res[1],res[2],res[3],rhon)
        1: fj = BGPS_FFORE_SIMPLE(R0,res[0],res[1],res[2],rhon)
        2: fj = BGPS_FFORE_MODEL2(R0,res[0],res[1],res[2],rhon)
        3: fj = BGPS_FFORE_SIMPLE2(R0,res[0],res[1],res[2],res[3],res[6],rhon)
        ELSE: message,"You've got a problem with mt.  RUM-TIME!!!"
     ENDCASE
     
     ;; DO vector
     R = d
     Z = d*0.d
     rho_hi = hi_density(R,Z)
     rho_h2 = h2_density(R,Z)
     rho_star = star_density(R,Z)
     
     CASE res[5] OF
        0: fj = BGPS_FFORE_MODEL(R0,res[0],res[1],res[2],res[3],rho)
        1: fj = BGPS_FFORE_SIMPLE(R0,res[0],res[1],res[2],rho)
        2: fj = BGPS_FFORE_MODEL2(R0,res[0],res[1],res[2],rho)
        3: fj = BGPS_FFORE_SIMPLE2(R0,res[0],res[1],res[2],res[3],res[6],rho)
        ELSE: message,"You've got a problem with mt.  RUM-TIME!!!"
     ENDCASE
     
     ;; Normalize the IN-PLANE radiation field, and place into
     ;; radfield structure
     rho /= rhon[0]
     radfield[j].rgal = rho
     rho_hi /= rhin[0]
     rho_h2 /= rh2n[0]
     rho_star /= rstn[0]
     
     oplot,R/1.d3,rho_hi,thick=1,color=cgColor('BLK5'),linestyle=2
     oplot,R/1.d3,rho_h2,thick=1,color=cgColor('BLK5'),linestyle=3
     
     oplot,R/1.d3,rho,color=colors[cj],thick=thick[cj],linestyle=0
     
  ENDFOR
  
  ;; oplot,R/1.d3,rho_star,thick=2,color=cgColor('BLK4')
  print,m4_stat(rho_hi)
  print,m4_stat(rho_h2)
  ;; print,m4_stat(rho_star)
  
  al_legend,/top,/right,textcolor=textcolor,$
            ['Case '+string(indgen(nmods)+1,format="(I2)")+$
             ' ('+cgSymbol('chi')+'!u2!dred!n = '+$
             string(chmin[0:2]/float(ndof[0:2]),format="(F0.2)")+')',$
             srho+'!dHI!n',srho+'!dH2!n'],$
            box=0,charsize=leg_cs,linestyle=[0,0,0,2,3],thick=[thick,1,1],$
            colors=[cnames,'BLK5','BLK5']
  
  ;; Build histogram of R_gal for the sources & overplot
  CASE rotc OF
     0: junk = KDIST(rb3.l, rb3.b, rb3.vlsr, RGAL=rgal, /NEAR)
     1: junk = KDIST_IAU(rb3.l, rb3.b, rb3.vlsr, RGAL=rgal, /NEAR)
     2: junk = KDIST_CLEM(rb3.l, rb3.b, rb3.vlsr, RGAL=rgal, /NEAR)
  ENDCASE      
  
  ;; Solar Circle marker
  plots,[R0/1.d3,10^(!y.crange[0])]
  plots,[R0/1.d3,2.2],/continue,linestyle=2,thick=5,$
        color=cgColor('Goldenrod')
  
  ;; Min R_gal marker
  lmin = min(rb3.l)
  print,'LMIN: ',lmin
  plots,[R0/1.d3*sin(lmin*!dtor),10^(!y.crange[0])]
  plots,[R0/1.d3*sin(lmin*!dtor),22],/continue,linestyle=4,thick=2
  
  yhist = float(HISTOGRAM(rgal/1.d3, min=0., max=20., binsize=0.2, $
                          loc=xhist, /L64))
  
  !y.type=0
  axis, yaxis=1,yrange=[0,149.99],yst=9,ytit='N sources per 0.2 kpc bin',$
        ytickformat="(I0)",charsize=1.2,/save
  oplot,xhist,yhist,thick=4,psym=10,color=cgColor('BLK4')
  axis,xaxis=0,xtickformat='blank_axis',xr=xrange,/xst;,xticks=3,xminor=5
  
  multiplot,/doxaxis
  



  bottompanel:
  ;;=========================================================================
  ;; BOTTOM PANEL -- SINGS COMPARISON
  
  restore,'ancillary/sings_spirals.sav',/ver
  ;; MAKE_RADIAL_PROFILE_8UM
  
  cgPlot,[0,1],[0,1],/nodata,yr=[0.1,150],/ylog,/yst,xr=xrange,/xst,$
       xtit='R!dgal!n [kpc]',ytickformat='exponent10',charsize=cs,$
       ;; ytit='Surface Brightness Normalized at R!dgal!n = 8.5 kpc',$
       ytit='Normalized Surface Brightness',$;xticks=3,xminor=5,$
       color=cgColor('Opposite'),xmargin=[8,2],ymargin=[4,1]

  ;; MWCOLOR = 'Green'
  ;; galcolors = ['Red','Orange','Opposite','Gray','Dodger Blue',$
  ;;              'Cyan',MWCOLOR]
  ;; gal_ls = replicate(0,7)
  
  MWCOLOR = 'Black'
  galcolors = [replicate('Black',3),replicate('BLK4',3),MWCOLOR]
  gal_ls = [2,3,0,2,3,0,0]
  
  FOR i=0L, n_elements(sings)-1 DO BEGIN
     
     index = fix(round(sings[i].t))-2
     color = galcolors[index]
     linestyle = gal_ls[index]
     print,color,linestyle
     
     cgOplot,sings[i].rad,sings[i].norm,color=color,thick=3,linestyle=linestyle
     
  ENDFOR
  
  restore,'./obsolete/irdc_dist_model/hyperion_radial.sav',/ver
  
  
  ;; pind = WHERE( pr GE 3.5 AND pr LE 8, np)
  pind = WHERE( pr GE 0 AND pr LE 15, np)
  oplot,pr[pind],int[pind],color=cgColor(MWCOLOR),thick=15
  
  al_legend,/top,/right,['T = '+string(indgen(6)+2,format="(I0)"),$
                         'R12 Model'],linestyle=gal_ls,charsize=leg_cs,$
            color=galcolors,delimiter='    ',$
            box=0,textcolor=textcolor,thick=[replicate(3,6),15];,linsize=0.6
    
  
  
  myps,/done;,/mp
  
END
